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Abstract. In two dimensions the Gross-Pitaevskii equation for a cold, dilute 
gas of bosons has an energy dependent coupling parameter describing particle 
interactions. We present numerical solutions of this equation for bosons in 
harmonic traps and show that the results can be quite sensitive to the precise 
form of the coupling parameter that is used. 
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1. Introduction 

The experimental realization of Bose-Einstein condensation (BEC) in dilute atomic 
gases has prompted a great deal of work on the theoretical description of such 
systems. Underlying many theoretical treatments is the notion that the condensate 
behaviour is governed, to first order, by the Gross-Pitaevskii equation (GPE). Recent 
experiments JjJ and proposals |2|, [| ^] for BEC in two dimensional (2D) traps have 
triggered interest in the properties of condensates in 2D. Although long wavelength 
fluctuations prohibit BEC in a 2D homogeneous system at any finite temperature g , 
the presence of a trapping potential alters the density of states sufficiently that finite 
temperature 2D Bose condensates may be created in trapped systems ||. Such 
condensates obey a modified form of the Gross-Pitaevskii equation, which contains the 
correct description of scattering in 2D. In a recent paper we derived this description 
based on an approximation of the many-body T-matrix in terms of the off-shell two- 
body T-matrix. In this paper we present numerical solutions of the 2D GPE for ground 
state and vortex states of a zero temperature gas of bosons in a harmonic trap, and 
contrast the results with other forms that have been suggested recently Q . We show 
that the detailed numerical predictions are quite sensitive to the precise form of the 
interaction strength used. 

In the following section we summarize the scattering theory needed to describe 
interactions in a 2D homogeneous BEC, before discussing in section ^| the various 
approximations by which these results may be applied to a trapped gas. Finally, in 
section ^ we present numerical results for each of these approximations and discuss 
the significant differences which can arise between them. 

2. The Gross-Pitaevskii equation and scattering 

The condensate wave function ip( r ) is given by the solution of the two dimensional 
Gross-Pitaevskii equation 

- |^ V V(r) + U trap (rWr) + ^o.9 2D (r) |^(r)| 2 ^(r) = ^(r), (1) 

where I4rap(r) is an external trapping potential (which is generally harmonic in 
present BEC experiments), No is the condensate population, and fj, is the chemical 
potential. The coupling parameter <?2d(i") appearing in the non-linear term describes 
the interactions between two condensate atoms. 

A collision between two atoms in momentum states \k) and |m) which produces 
a transition to states \i) and |j) is described by the T-matrix element (ij\T(E)\km) , 
where E is the energy of the collision. The T-matrix is obtained as the solution of a 
Lippmann-Schwinger equation or equivalently via a summation of ladder diagrams 0. 
In three-dimensional (3D) systems the coupling parameter in the GPE is often taken to 
be the zero-energy, zero-momentum limit of the two-body T-matrix, which describes 
the scattering of particles in a vacuum. This gives (00|T 2 b(0)|00) = (73D = ^h 2 azT>/m 
where CI3D is the s-wave scattering length ||. This is, however, merely an approximate 
description since the scattering of two condensate particles actually occurs in a 
medium consisting of the surrounding particles rather than in a vacuum. Instead 
the collision is properly described by a many-body T-matrix element (00|Tmb(0)|00) 
which incorporates the effects of the surrounding atoms on the scattering process. In 
3D the many-body T-matrix leads to a relatively small correction to the two-body 
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T-matrix approximation (of relative order (na^) 1 / 2 at T = 0) and for many purposes 
it is sufficient to neglect many-body effects in the GPE. In two dimensions (and 
lower), however, the two-body T-matrix vanishes in the zero-energy, zero-momentum 
limit [ flCfl , and many-body effects are therefore of much greater importance and 
contribute even at leading order. 

In a recent paper we demonstrated how the many-body T-matrix can be 
approximated by the two-body T-matrix evaluated off the energy shell. The coupling 
constant which appears in the GPE in a homogeneous 2D system was shown to be 

<?2D = (00|T MB (0)|00) = (00|T 2b (- M )|Q0>. (2) 

Using the expression for the off-shell two-body T-matrix found in reference Jl0| , gives 
the following form of the coupling parameter in a homogeneous system 

Airh 2 1 

In (fima 2 ,^ / Ah 2 

where a2D is a two-dimensional scattering length. For a 2D gas of hard spheres of 
radius a we have <22D = ae 7EM , where 7em ~ 0.577 is the Euler-Mascheroni constant 
which we have absorbed into the definition of a2D here for convenience^. In practice, 
a 2D gas is created by trapping atoms very tightly in one dimension (the z axis) such 
that the motion in the z direction is effectively frozen out. The effective 2D scattering 
length for such a gas is given by Jlf|, [T^] 

«2D = 4 A \^-l z exp ( -Vtt— ] , (4) 
V B V a 3D J 



^2D _ c. 2 /Ati2\> W 



where B sa 0.915, a3D is the 3D s-wave scattering length, and l z = ^Jh/2mLd z is 
the typical width of the system in the z direction. The 2D scattering length therefore 
depends not only on the 3D scattering length, but also upon the degree of confinement 
in the z direction. 

Equation (||) shows that the coupling parameter for a 2D homogeneous Bose gas 
depends on the chemical potential of the system (and hence the density) as well as 
the 2D scattering length. This is in contrast with the case in 3D where the coupling 
parameter depends only upon the scattering length to first order. 



3. The GPE in trapped 2D systems 

The expression for the 2D coupling parameter shown in equation (||) was derived for 
a homogeneous system, and the correct application of these results to the case of a 
trapped gas is the main objective of this paper. In a previous paper we provided 
solutions of the 2D GPE in a trap using the homogeneous coupling parameter of 
equation (||) in order to illustrate the effect of the energy dependence of g^D ■ In this 
paper we focus on a more accurate description of the scattering in inhomogeneous 
systems. 

We consider a 3D Bose gas confined tightly in one dimension and weakly in the 
remaining two dimensions on a length scale ftrap- A collision between two condensate 
particles will typically occur over some length scale £ C oll- Provided that t co \\ is much 
smaller than ftrap we can use a local density approximation. 

| Note that this definition of ci2D differs sightly from that in our earlier work such that our a2D 
here equalsa2De 7EM in the notation of reference Jip| j . The definition used here simplifies the form of 
equation (^). 
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We can introduce the length scale £ co \\ by the following simple argument. We 
model the pair wave function of two atoms in the medium by that of a single particle 
with the reduced mass moving in a potential which consists of a circularly symmetric 
box of radius L and a hard sphere of radius a located in the centre of the box. For 
s-wave scattering the wave function is solved by 

i>(r) =A J (kr)+B N (kr), (5) 

where Jo and No are Bessel functions of the first and second kind respectively. In 
the zero-energy, zero-momentum limit we get ip( r ) = ^-o + £?ohi(r). Applying the 
boundary conditions that the wave function vanishes on the radius r = a and reaches 
the asymptotic value \ at the edge of the box gives 

iKr) = xPrll\ for a < r < L. (6) 
m(L/a) 

The extra energy caused by the curvature of this wave function resulting from the 
presence of the scattering potential is 

a* * /WjfA-iQs!!* (7) 



2m J a 2m ln(L/a) 

This energy depends upon the size of the box L, which is indeed the length scale 
relevant for the scattering of two particles in 2D. The scattering of two particles in 
a many-body system should obviously not depend on the size of the system as a 
whole when L becomes large, and so we must interpret L as the physically relevant 
length scale £ C oii- The appropriate length scale over which a many-body wave function 
changes is the healing length £h, given in homogeneous Bose condensed systems by 
&h = ft/ V2m52D"-o = h/y/2mfi, and so it is this which must be used in equation ([?]). 
Since A^ 1 X 1 2 corresponds to the condensate density n , this leads in the homogeneous 
limit to a pair interaction strength of the form of equation (||) . 

The same argument can be applied straightforwardly to trapped gases if the 
condensate density varies slowly on the scale of the healing length, which is true 
except in the surface region. In this case, many-body effects cause the pair wave 
function to reach its asymptotic value on a length scale equal to the local healing 
length lh = %/ y / 2mg2i){^)no(r). The two-body interaction strength is therefore given 
from equation (||) by replacing fi with no(r)g2D(r) producing a density-dependent 
effective interaction. Such density dependent coupling parameters are expected from 
the results of density functional theory |Tj| which predict that the energy of the system 
is a functional of the density only. 

In an inhomogeneous system the density is spatially dependent and thus the 
coupling parameter is also spatially dependent. In terms of the condensate wave 
function the density is given by no(r) = No\i^(t)\ 2 . Equation (0) now gives for the 
coupling parameter the result 

4irh 2 1 
32D(r) ~ m ln(iV |^(r)|2 52D (r)ma2 D /4^)- W 

An approximate solution to this equation may be found by iteration, giving 

Wr) = [WNoWWatt] 1 + O ( . (9) 

m \ \n(n a^ D ) J 

The first order term in this expansion agrees with the form of coupling parameter 
proposed by Kolomeisky et al. @, IT3] who used the renormalization group to analyse 
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a 2D homogeneous Bose gas. Earlier work by Shevchenko [[15] also proposed such 
a coupling parameter based on the work of Schick ]l6[ |. More recently Tanatar et 
al. Jl7t have also made use of this approximation. Unfortunately the expansion in 
equation (Q) is not valid for realistic systems since the higher order terms are not in 
general negligible, as will be shown in the following section. 

A more accurate procedure would be to solve equation (||) numerically for 
<72D(«o(r)) and use this exact solution in solving the GPE. In the following section 
we present results which suggest that this accurate solution may be necessary in some 
circumstances. 



4. Solutions to the Gross-Pitaevskii Equation 

In this section we present numerical solutions of the 2D GPE for a Bose condensate 
trapped in a circularly symmetric, harmonic potential characterized by an angular 
frequency uj. The various approximation schemes for the coupling parameter discussed 
in the previous section will be compared to the more accurate <72d(i") found from 
the numerical solution of equation (||). In order to illustrate our results, we choose 
as our parameters to = 2tt x 100Hz and a2D = 6nm. This 2D scattering length is 
approximately the 3D scattering length found for 87 Rb, and so from equation (jl]) this 
corresponds to a situation where l z ss a^, and hence the two-dimensional nature of 
the scattering is important . For these parameters the results we obtain correspond 
to healing lengths between O.l^trap and 0.5£trap (for the very low /i solutions) at the 
centre of the trap. These healing lengths are sufficiently small for the local density 
approximation to be valid. 

We will solve the GPE for three different approximations for the coupling 
parameter. The simplest case restricts g2D(r) to be constant everywhere and 
determined by equation (|^) as in our earlier work Spatial variations can most 
simply be introduced by using the first term of the expansion in equation (^), which 
corresponds to the results of Kolomeisky et al. [|[ [lij. Finally, the most accurate 
approximation is obtained using the full numerical solution of equation (||). Figure [l] 
provides sample solutions of the ground state wave functions and coupling parameters 
calculated using these three different approximations for the same chemical potential. 
It can be seen that, although the wave functions are fairly similar, the coupling 
parameters behave quite differently. The results for the constant coupling parameter 
agree well with the predictions of equation (||) , and differ significantly only towards 
the edge of the condensate. Of course, this is to be expected since fj, ~ n.off2D(r) 
near the centre of the trap. The energy contribution due to interactions for atoms on 
the edge of the condensate is greater with the constant coupling parameter than with 
either of the spatially varying parameters, and hence the constant parameter wave 
function has a lower amplitude in these surface regions. 

Figure H does show a large difference between the two spatially dependent coupling 
parameters, however, especially in the central region where the condensate density is 
greatest. The coupling parameter of equation (^) is greater by about a third at the 
centre of the condensate than the full expression of equation (@), and remains larger 
throughout. This arises because the expansion of equation (Q) does not converge 
sufficiently rapidly. Indeed, for the case illustrated here, the second order term in the 
expansion (which is negative) reaches a magnitude of approximately —1 at the centre 
of the trap (in the units used in figure |l]b) . 

Such large differences in the coupling parameters can lead to problems when 
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Figure 1. (a) Wave functions V( r )> an d (b) coupling parameters g(r) calculated 
for fi = 25frur. The dash-dot lines correspond to equation ([?]), the dashed lines 
correspond to equation (9h, while the solid lines represent the full numerical 
solution of equation (fcj). 
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calculating related quantities, such as condensate populations. Figure ^| shows the 
results for condensate populations versus chemical potentials calculated with each 
of the coupling parameters. As can be seen the predictions obtained using the 
approximation of equation (^) underestimate the condensate numbers by roughly 20 
per cent compared to the full numerical solution of equation (||) . This is a consequence 
of the greater strength of the coupling parameter which occurs in this approximation. 
Agreement between equations (^) and (Q) is poor even though the usual criterion for 
a dilute gas (nodfrj <C 1) is obeyed (noa^D ^ s °f the order of 10 for the situation 
illustrated in figure [l]). Indeed, in order to apply the approximation in equation (^J) 
we require that ln(noa2 D ) 3> 1 (while ^o a 2D 1)> which is a much more stringent 
criterion, and one that is experimentally unfeasible requiring at least no a 2D 10 -20 . 
For this reason the use of the full expression in equation (||) is necessary to simulate 
real 2D gases. 

The spatially-constant approximation to the coupling parameter gives compara- 
tively much better agreement with the full numerical solution in figure although it 
also underestimates the condensate number by about 5 per cent. It would appear from 
these results therefore that if a simple analytical approximation of <72D is required to 
estimate bulk properties of the system then the spatially constant approximation is 
preferable to equation (|^) . Of course the spatially constant approximation does not 
deal with the boundary regions of the condensate well, and so for properties dominated 
by edge effects the better analytical approximation is likely to be that of equation (|9|) . 

The wave functions for vortex states can also be obtained in 2D, if we assume a 
solution of the form 

^(r) = <f>(r)e iKe , (10) 

where 8 is the angle around the vortex core, and n is an integer. The phase of ip 
therefore wraps around by 2itk as the range of is traversed. The energy per particle 
(in a non-rotating frame) for a condensate with wave function ip is given by the 
functional 

.(11) 



E[ip] = J dr 



|^|W(r)| 2 + Wr) Wr)| 2 + ^f^l^(r)| 4 



Creation of a vortex in the centre of the trap comes at the cost of increasing the 
contributions from both the kinetic energy and the trapping potential terms in the 
energy functional, although the interaction term is reduced by virtue of a lower central 
density. The single vortex state can be made energetically favourable by rotating 
the trap at a frequency O such that E[ip K =o] becomes less than E[tp K —i] — ilNofon. 
The point at which this occurs is known as the thermodynamic critical frequency 
Q c , and this is shown in figure I for the various forms of #2D- The 2D critical 
frequency is substantially lower than for a 3D gas with the same scattering length 
(<J3D = G2d) due principally to the much higher interaction strength which occurs in 
2D 0. In figure ||b it can be seen that the effect of a density-dependent coupling 
parameter is to reduce the critical frequency as compared to the constant parameter 
case. This is to be expected since the appearance of a vortex lowers the mean density 
of the condensate, which decreases the coupling parameter calculated from either 
of equations or (^|). The lower coupling parameter means a greater saving in 
the interaction energy term when a vortex is created and therefore decreases the 
critical frequency. The saving in the interaction energy is due principally to the 
reduction of the density in the centre of the condensate where the density (and 
hence coupling parameter) is greatest in the ground state. Because the ground 
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Figure 2. a) Condensate numbers as a function of chemical potential in 2D 
for various coupling parameters, b) The percentage differences in the condensate 
populations as compared to the full numerical solution. The line styles correspond 
to those used in figure |l . 
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Figure 3. a) Critical frequency versus condensate number for 2D and 3D 
condensates. The upper curve corresponds to the 3D case with = ci2D> 

whilst the lower curves represent the various 2D results, b) A detailed view of 
the critical frequencies in 2D using the three forms of the xoupling parameter 
discussed. The line styles correspond to those used in figure 111 

state coupling parameter is much larger in the approximation of equation (|]), the 
interaction energy saving is also greater. The critical frequency is therefore lower in 
this approximation compared to the results obtained using the coupling parameter of 
equation (||). In contrast to the spatially dependent 52D cases, the spatially-constant 
coupling parameter calculated from equation (|^) increases when a vortex is formed, 
due to the greater chemical potential of the vortex state, and so the critical frequency 
is higher in this approximation. 
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5. Conclusions 

In this paper we have applied previous results obtained for the many-body T- 
matrix in a homogeneous condensate to the more currently relevant problem of a 
trapped condensate, by means of a local density approximation. This leads to a 
spatially dependent coupling parameter appearing in the non-linear term of the Gross- 
Pitaevskii equation. We have shown that results obtained using the full numerical 
solution to the coupling parameter can differ substantially from the simple first 
approximation obtained via a series expansion. The form of the coupling parameter 
in this approximation is the same as that presented in recent work by Kolomeisky et 
al. (|, , Tanatar et al. [jl7| , and Shevchenko jL5| , and is closely related to the work 
of Schick flq| . However, this approximation is only valid in the limit ln(noa2D) 1 
(while /loa^D < -0 which is not experimentally relevant. Our results indicate that the 
full expression of equation (||) may be needed to model current experiments accurately. 
Corrections to equation (||) are expected to be of order (na^)/ ln(na| D ), and the limit 
where this parameter is small should be experimentally relevant. 

Agreement with the full numerical solution for the coupling parameter is found 
to be substantially better if it is approximated using a spatially constant (but energy 
dependent) parameter as in the homogeneous limit. It would seem from the results 
presented here that if an approximate analytical form of the coupling parameter is 
required (for deriving approximate Thomas-Fermi wave functions for example) then 
the spatially constant form of <72D given in equation (||) is preferable to the expression 
of equation (^J) for many purposes. 
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